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The package CosmoLib is a combination of a cosmological Boltzmann code and a simula- 
tion toolkit to forecast the constraints on cosmological parameters from future observations. 
In this paper we describe the released linear-order part of the package. We discuss the 
stability and performance of the Boltzmann code. This is written in Newtonian gauge and 
including dark energy perturbations. In CosmoLib the integrator that computes the CMB 
angular power spectrum is optimized for a £-by-£ brute-force integration, which is useful for 
studying inflationary models predicting sharp features in the primordial power spectrum of 
metric fluctuations. As an application, CosmoLib is used to study the axion monodromy 
inflation model that predicts cosine oscillations in the primordial power spectrum. In con- 
trast to the previous studies by Aich et al and Mccrburg et al, we found no detection or 
hint of the osicllations. We pointed out that the CAMB code modified by Aich et al does 
not have sufficient numerical accuracy. CosmoLib and its documentation are available at 
http : //www. cita.utoronto . ca/~zqlmang/CosmoLib 
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INTRODUCTION 



The hot big bang model and the cosmological perturbation theory, where the physical metric is 
perturbed around the spatially homogeneous and isotropic Friedmann-Robertson- Walker (FRW) 



metric 



have led to a remarkable success in interpreting the plethora of observational data of 
the last two decades |5l— tldt] . Observations of the temperature anisotropy in the Cosmic Microwave 
Background (CMB) have been playing an essential role in building the standard cosmological model 



and measuring its parameters 



10 



111 ]. In order to maximize the usage of the observational data, 



one would like to compute the theoretical prediction on the CMB anisotropy for a given model as 
accurately as possible, wit : 



years such as CMBFAST 



i__tolerable time consumption. Computation tools developed over the 
Q, CMBEASyQand CLASS Q, Q are capable of 



12j, CAMB 

computing a CMB angular power spectrum to percent-level accuracy within a few seconds on a 
modern desktop personal computer. 

The crucial technique used in all the fast CMB codes to date is the line-of-sight integration 
approach ,[li| and an assumption that the primordial power spectrum of metric perturbations 

is smooth. (Here for readability we focus on the comoving curvature perturbations and temperature 
anisotropies, although the same arguments can be as well applied to the tensor perturbations and 
CMB polarization.) A CMB code first computes the radiation transfer function by solving the 
linear-order Boltzmann equations and using the line-of-sight integration method, then convolves 
|A^| 2 with the primordial power spectrum V(k) to obtain the CMB angular power spectrum Cg. 
The smoothness assumption allow us to compute only a few tens of multipoles spanning from 
^min = 2 to £ max ~ a few thousands and interpolate the remaining C/s. Furthermore, since V(k) 
is assumed to a smooth function, sparse sampling of the radiation transfer function has been 
implemented for the integration of each Cg. 

A smooth primordial power spectrum is a prediction of the simplest single-field slow-roll infla- 



tion models 



19M23I] . However, local signatures in the primordial power spectrum that makes it 



deviate from smoothness can arise in various alternative models, for instance, when the inflaton 
potential has sharp features 



ton evolution 
inflation [3d . 



J3 



24, 



251 ]. when there is a transition between different stages in the infla- 



271 ] , when more than one field is present [28|, |29( , from particle production during 



311 ] . modulated preheating 



odromy in the extra dimensions 




, or, more recently, in models motivated by mon- 



341 ] (see also [35J). These features represent an important window 



on new physics because they are often related to UV scale phenomena inaccessible to experiments 
in the laboratory. For these models, the CMB angular power spectrum is not necessarily smooth, 
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and therefore needs to be computed at each multipole without interpolation. This increases the 
computing time by a factor of a few tens. Moreover, for the numerical- integration of each Ce, the 
sampling frequency in the wavenumber k often needs to be increased, again, by a factor of a few 
tens. The required sampling frequency in k is model-dependent. It is determined by the larger 
between the minimum width of the features in the primordial power spectrum and the minimum 
width of the oscillations in the radiation transfer function. 

To keep track of the features in the primordial power spectrum, one can modify standard 
CMB codes by naively doing an l-by-l brute-force calculation with increased integration sampling 
frequency in k. However, in the case where the features in the primordial power spectrum are really 
sharp (Sink < 0.01), this naive modification increases the computing time by a factor of ~ 10 3 (a 
few tens in i sampling times a few tens in k sampling). Moreover, the memory that is required to 
store all the transfer functions and tables of spherical Bessel functions can be too large for most 
desktop personal computers. One of the purposes of this paper is to introduce a more optimized 
algorithm to treat these problems. In fact, apart from increasing the sampling frequency, that 
cannot be avoided, all the other problems can be significantly alleviated by using the recurrence 
relation of spherical Bessel functions. An optimized algorithm, which we detail in Section IIIIl is 
< 10 2 times slower than the standard algorithm for the smooth-'P(A;) case. This new algorithm 
has been implemented in the CosmoLib package, a self-contained package that we developed to 
compute cosmological perturbations, CMB angular power spectra, and the forecast constraints on 
cosmological parameters from future cosmological surveys using Fisher matrix analysis and Monte 
Carlo Markov Chain (MCMC) calculation. In particular, the cosmological surveys that we consider 
are CMB, large scale structure (LSS) and supernovae (SN). 

In addition to the enhanced CMB integrator, CosmoLib has a few other features that are 
complementary to the other publicly available Boltzmann/CMB/MCMC codes. For instance, the 
MCMC engine in CosmoLib has a modified rejection rule that allows the proposal density (the 
probability of random-walking to a new point in the parameter space) to periodically depend 
some parameters. This is useful when one considers a likelihood that depends on some per iodic 
parameter. This happens, for instance, in the context of inflation from axion monodromy 



3, 3- 



381 ] . where the oscillations in the predicted power spectrum depend on a free phase. Moreover, 



CosmoLib treats the dark energy equation of state (EOS) w(a) and the primordial scalar and 
tensor power spectra V s (k) and Vt(k), as free functions, which can be either chosen from a list 
of build-in models or defined by the user. This makes CosmoLib a convenient tool to study non- 
standard parametrizations of dark energy EOS and primordial power spectra. Finally, CosmoLib is 
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TABLE I. Comparison between CMB Codes 





tAMB 


CLAbb 


L MBbAb Y 


CMB quick 


CosmoLib _ 


Language 


b 90 


c 


C++ 


"ft T ^ i 1 _ j * „ 

Matnematica 




gauge _ 


syn. 


syn. /Newt. ^ 


syn. /gauge- inv. 


Newt. 


Newt. 


open/close universe 


Yes 


AT 

INo 


INO 


No 


No 


massive neutrinos 


Yes 


Yes 


Yes 


Yes 


AT 

No 


tensor perturb. 


Yes 


Yes 


Yes 


Yes 


"Win 

Yes 


{' ■ I '1 TI T' . 1 

CUM isocurvature mode 


Yes 


Yes 


Yes 


Yes 


Yes 


dark energy perturb. 


Yes 


Yes 


Yes 


INO 


Yes 


2 

nonzero c s h 


Yes 


Yes 


Yes 


No 


Yes 


dark energy EOS. 


constant 


w + w a (l - a) 


arbitrary 


-1 


arbitrary 


non-smooth primordial power 


No 


No 


No 


No 


Yes 


MCMC driver 


Yes 


No 


Yes 


No 


Yes 


periodic proposal density 


No 


No 


No 


No 


Yes 


data simulation 


No 


No 


No 


No 


Yes 


second-order perturb. _ 


No 


No 


No 


Yes 


Nof 



a Here we do not include CMBFast, which is no longer supported by its authors or available for download. 

This refers to CosmoLib Version 0.2. 
c CosmoLib is a mixture of Fortran and C codes. The main part is written in Fortran. 
d syn.: synchronous gauge; Newt.: Newtonian gauge; gauge-inv.: gauge- invariant variables. 
e Newtonian gauge is implemented in CLASS version 1.3. 

A second-order perturbation code is used to study the CMB non-Gaussianity. 
6 The second-order part of CosmoLib is not released with this paper. 



written in Newtonian gauge (also called Poisson gauge) |39r44l[]. while many other codes are mainly 
developed in synchronous gauge (see e.g. [^j]). This is a plus-and-minus point. We found that our 
Newtonian-gauge Boltzmann code is slightly slower than the codes written in synchronous gauge. 

Newtonian gauge. For 
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551 ] . Implementing 



However, many theoretical works in the literature have derived equations in 
instance, second-order Boltzmann equations have been derived in this gauge 
these equations in a code already in Newtonian gauge would be much easier. To conclude the 
discussion, we list the differences between CosmoLib and other publicly available CMB codes in 
Table E 

As an application, CosmoLib is used to study the "hints" of cosine osicllations in the primor- 
dial power spectrum that was recently found in Refs. jsfl In an accompanying paper 58]. 
CosmoLib is applied to forecast the constraining power of future CMB and galaxy survey data on 
the primordial power spectrum from inflation, with an emphasis on models generating features in 
the power spectrum. 
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This paper is organized as follows. In Section [II] we introduce the Boltzmann code in Newtonian 
gauge and discuss its stability and performance. Section IIIII details the algorithm used in the 
enhanced CMB integrator. In Section IIVI we introduce the forecast technique and parameter 
estimation methods. Section IVl concludes. 

Throughout this paper, unless otherwise specified, repeated indices are summed over. Greek 
indices run from to 3. Latin indices run from 1 to 3, that is only over spatial dimensions. We 
use natural units c = h = 1 and the reduced Planck Mass M p = 1/\/8ttGn = 2.43 x 10 18 GeV. 



II. COSMOLIB IN NEWTONIAN GAUGE 



A. The Background Solutions 



Let us start discussing the background solutions. We consider a flat FRW metric ds 2 = 
a 2 (r)(— dr 2 + dx l dx l ), where a is the scale factor and r is the conformal time. The normaliza- 
tion of a is arbitrary. We normalize it such that a = 1 today. CosmoLib uses the e-fold number 
N = In a as the time variable. The physical Hubble expansion rate is defined as H = da (f T . Its 
present value is denoted by H = 100/ikms" 1 Mpc -1 . 

We assume a universe with cold dark matter (labeled with a subscript c), dark energy (labeled 
with a subscript A), baryons (labeled with a subscriber b), radiation (labeled with a subscript 
7), and 3 species of massless neutrinos (labeled with a subscript v). For a component X (X = 
b, c, 7, v, A) the background density is denoted as px, and the background pressure px- The present- 
day fractional energy density is written as £lxo ■ Dark energy is assumed to be a perfect fluid with 
known equation of state (pressure to density ratio) w(a) and a constant sound speed in its 
rest frame. The users can either choose w(a) from a list of build-in models or define their own 
w(a) functions. The build-in models of w(a) include the cosmological constant model w(a) = — 1 



59( | . a constant EOS 10(a) = wq, a linear function w(a) = wq + w a (l — a) [6CJ, [6l|, and a general 



three-parameter parametrization for the minimally coupled quintessence/phantom models 



62J 



For a given w(a) the background solutions are 

,AT 



a 



e 

3H 2 M 2 n c0 a' 



, Pc 



0, 



Pb = 3H£Mft b0 a- 6 , p b = , 
p 7 = 3H$Mfa 70 a~ 4 

Pv 



P 1 



1 

3ifo M p^^oa" 4 , Pv = ^Pu , 



(1) 



PA = 3HQMptt A0 a 3 exp 
H 



o 



A 



w(a)dN 



p\ = w(a)pA , 



1 Pc + Pb + P 7 + Pu + PA 



We will also use the derived quantities Qx( a ) = Px I '(3if 2 M 2 ) (X = 6,c,7, i/, A), i? = (3p&)/(4p 7 ), 
and 



din if 

diV 



1 + 



PA + P 7 + Pi, 



+ Pb + PA + P 7 + Pv 

N 



The conformal time r can be related to the scale factor a = 

" a da 



T 







if a 2 ' 



(2) 



(3) 



The electron number density n e (a) is obtained using RecFast version 1.5 |63j, [64|, which has 
been incorporated into CosmoLib. We denote the differential optical depth (increment of optical 
depth per dN) as 



dn 



K N 



n e a T 



(4) 



dN H 

where gt = 6.652 x 10 _25 cm 2 is the Thomson scattering cross section. The baryon sound speed 



c 2 b (a) is obtained by solving the differential equations (68-69) in Ref. 



With these background solutions in hand, now we can write down the governing equations for 
scalar perturbations. 



B. Scalar Perturbations 

The metric in the (generalized) Newtonian gauge can be written as 

ds 2 = a 2 (r) {-(1 + 2$)dr 2 + Uidx'dr + [(1 - 2f)<% + h {j ] dx l dx j ] , (5) 

where diOJi = 0, ha = and dihij = 0. The vector perturbation cjj decays in an expanding universe 
and hence it is set to zero in CosmoLib. The tensor perturbation hij is gauge-invariant and its 
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governing equations are identical in all gauges. Thus, we will only focus on the scalar perturbation 
equations that in CosmoLib differ from those in many other Boltzmann codes. 

The linear-order relative density perturbation of X is denoted by 5x = Spx/px, and the linear- 
order velocity vx- Unless otherwise specified, 5x and vx are all defined in Fourier space, that are 
functions of r and the wave vector k. 

The radiation relative temperature fluctuation AT/T from direction n seen by an observer at 
position x is expanded as 



AT , f d 3 k 



/7ol <sv ^ j 

(2^3 E^E 2 q t(^ ™)H) y 4(2zti) yr(n)eikx ' (6) 

where Y™ are the spherical harmonic functions. Note that the moments 7 (£, m) are functions 
of the wavenumber k and the conformal time r. The energy density fluctuation and velocity of 
photons are related to the moments I = and £ = 1 via 

<5 7 = G 7 (0,0); u 7 = ^e 7 (l,0). (7) 

The neutrino moments Q v (£,m) are defined in the same way, by replacing the subscript 7 with v. 

For the polarization of radiation, the Stokes parameters Q, U are expanded using the spin-2 
harmonics ±2^ m 13] 

d 3 k 

£=0 m=-2 

where E(£,m) and B(£,m) are functions of the wave vector k and conformal time. 

The linear-order Fourier modes are decoupled. The Fourier-space variables to be evolved are 
V N = d9/dN, 5 b , v b , 5 C , v c , 5 A , 9a = (l + w)v A , 9 7 (£,0) (£ = 0, 1, 2, .., £ max , 7 ), e„(*,0) (£ = 0, 1, 
2, .., £ maXiU ), E(£,0) {£ = 2, 

•^maxji?)- The truncations ^max^j 4iax,i/ and ^ max ,E are adjustable 
integers. In CosmoLib their default values are taken to be 14, 12, 14, respectively. Without loss of 
generality we choose the azimuthal direction (the z-axis direction that is used to define Y^ m (n)) 
to be parallel to k. 

The gravitational potential «3> can be obtained from the Einstein equations 



(Q±iC/)(x,n,r) = J ^ £ £ [E(£, m) ± B{£, m)] H)^ ^ + 1} [±2*7" (n)] e* k ' x , (8) 
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40. 



65] 



* = * " [^7 Q 7( 2 > °) + n„e„(2, 0)] , (9) 



where we have introduced the reduced wavenumber 



k " s J? ■ < 10 > 



8 



Note that kn, ^ 7 , are functions of time. We do not treat $ as an independent- variable. Instead 
we view it as a function of the variables \&, O 7 (2,0) and 0„(2,O). 

The close set of first-order differential equations including all the truncation schemes is: 



dN 
d5 c 

dN 
dv c 
dN 
d5 b 

dN 
dvb 
dN 
d5 A 
dN 



d9A 
dN 



d9 7 (0,0) 
dN 

de 7 (i,o) 

dN 

d6 7 (2,0) 
dN 

dG 7 (£, 0) 
dN 



de 



dN 

rfe„(o,o) 

dN 

de v (i,o) 

dN 
d&„(£,0) 



dN 



de 



7V < -max,n 



dN 

dE (2,0) 
dN 

dE{£, 0) 
dN 



-k H v c + 3*at , 

-V e + fcff* , 

-k H v h + 3^ N , 

-v b + k H {$ + c 2 S)b S b ) - ^ 



^ 6 --e 7 (i,o) 



-3 (c 



W 



)5 A -9 



w 



dw/dN 



^-k H 6 A + 3(l + w)V N , 



w + c 



s,A 



w 



dw/dN 
3(l + w) 



3(1 + w 

0a + k H [c 2 sA S A + (1 + w)S] 



1 

~3 
k H 

k H 
k H 



k H @ 7 (i,o) + m N , 



e 7 (o,o) - -e 7 (2,o) + 4$ 



3 J \ 3 



|e 7 (i,o)-^e 7 (3,o) 



K N 



l 07 (2,O) + ^(2,O) 



(11) 

(12) 

(13) 
(14) 

(15) 

(16) 
(17) 

(18) 
(19) 

(20) 



£ + 1 



2£ 



_e 7 (£_l )0 )___G 7 (£ + l,0) 



k n Q^(£,0) (2<£<£ max , 7 ) , 



7V c -max,7; 0) 2^ maXj7 + 1 



2f 



max, 7 

^ H e„(i,o) + 4f w , 



+ 1 



-fcf/9 7 (4nax, 7 - I) ) - ( K JVH ^jr: j @ 7 (^max,7, 0) 



k H 



k H 



21-1 

0) 2£ maXjJy + 1 

~~ ~~ 2f 



6„(0,0) - -6^(2,0) +4$ 
5 

e^i-1,0)- -±Ie^ + i,o) 



k H @ 



v V l max,f 



1^£(3,0) - k n 



k H 



7 

-^1,0,2 

2£-l 



2£ + 3 
-1,0)- 

2 



+ 1 



aHr 



(2 £ <^ fniax,!/) > 



6f(^max,i/! 0) , 



5 £(2,0) + ^6 7 (2,0) 



^-1,0) -^^ + 1,0) 



k n E(£,0) (2<£< 



(21) 
(22) 

(23) 
(24) 

(25) 

(26) 

(27) 
(28) 



9 



+ 1 



dN 



2C 



max,E 



\k H E(£^ E -1,0)- U N + Ca a X ^ r +1 ) #(W,£, 0) , (29) 

dw/dN \ 6a_ 
~ W + 3(1 + w) J fa 



(3<a-1)*a + 9U 



fcu , . 3 / dOv(2,0) „ de„(2,0) 



(30) 



In the radiation and neutrino hierarchy equations (| 18ti29j) we have used the Clebsch-Gordan 



coefficients Ki ma , which are defined as 



17] 



— ^ , it i > max{|m|, \s\, 1) ; 

, otherwise. 



(31) 



These equations are already written in the form that can be directly implemented into a generic 
first-order ordinary-differential-equation (ODE) solver. For a derivation of these equations, see 



Refs. 
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65, 



661 ] . (The change of time variable from r to N can be done straightforwardly 



using d/dr = aHd/dN .) The initial conditions can be found in Ref. 



approximation we follow Ref. 



401 ] . For the tight-coupling 



671 ] . where the obvious typos in Eqs. (15-16) has been fixed. Since 



these treatments are identical to the original source, we do not repeat the discussion here. The 
interested readers are referred to these references for the governing equations and their derivation. 

CosmoLib allows the user-input w(a) to be a phantom-crossing function, that is a function 
crossing the line w = —1. In this case we force d5\/dN and d9\/dN to be zero around the 
proximity of the phantom crossing. This is an approximation. Exact treatment requires input of 



at least one more degree of freedom 
Ref. 



6814701] . which cannot be implemented in a generic code. In 



661 ] the reader can find an alternative treatment that works better for multiple scalar field 
models. 

Equation (130p is the key equation that guarantees the numerical stability of the code (even for 
isocurvature initial conditions). It is obtained by subtracting the ii components of the perturbed 
Einstein equations (pressure perturbations) from the 00 component (density perturbations). This 
particular combination of the Einstein equations has been applied in the numerical code CMBquick 



50, 



711 ]. which assumes that dark energy is a cosmological constant and ignores the baryon sound 



speed. Eq. (|30[) is a generalized version that includes dark energy perturbations and a nonzero 
baryon sound speed. 

We can use the energy constraint (00 component of the perturbed Einstein equations, that is 
^G*oo = STqq) and the momentum constraint (Oi-component of the perturbed Einstein equations, 
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OT, 



|ST 00 | 
5G ool 
|5T 0i l 
|5T 0| - 5G„J 



k = 0.1 h Mpc" , adiabatic 



0.1 
0.01 
0.001 
0.0001 
1e-05 
1e-06 
1e-07 
1e-08 - 
1e-09 

-30 



|5T 00 | 
|5T 00 - 5G D0 | 
|5T 0i l 
|5T 0| - 5G„J 



k = 0.0001 h Mpc , adiabatic 




1ST 



[ST 00 | 

|5T 0i l 
iST 0i - 8G„,| 



k= 0.01 h Mpc" 1 , adiabatic 




FIG. 1. Testing the energy constraint (00-component of the perturbed Einstein equation) and momentum 
constraint (Oi-component of the perturbed Einstein equation). The cosmological parameters used here are 
Qbah 2 = 0.022, ftcoh 2 ~ 0.1128, h = 0.72. In the lower-right panel the CDM-isocurvature initial conditions 
are used, while in the other panels we have used adiabatic initial conditions. 

that is SGoi = 5Toi) to estimate the numerical error of the code. As shown in shown in Figure [H 
the relative errors are < 10 -4 for a wide range of scales and different initial conditions. 



III. CMB ANGULAR POWER SPECTRA 
A. Algorithm 

Optionally CosmoLib can compute the CMB angular power spectrum for each multipole t by 
brute force, i.e., without interpolation. The angular spectrum for the temperature anisotropies is 
given by 

C t = J \A k e \ 2 V s (k)dlnk , (32) 
where is the temperature transfer function given by the line-of-sight integration 

[■TO 

A k £ = / S(Kr)u[k(T -T)}dr , (33) 
J o 
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Li, i 
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/ = 300 



5 6 
In (k/Ho) 



FIG. 2. The temperature transfer function A{? for a fixed I = 300. A typical sampling scheme is shown by 
the red solid triangles in the upper-right panel, which zooms-in part of the figure. 



where je is the spherical Bessel function and tq is r at redshift zero. The source S(k, r) can be 



12|,[l7]. 

3, 



computed from the perturbations ty, <£, Sx, vx (X = c, b, A, 7, u), G 7 (2, 0) and 0^(2,0) 
In Ref. jl?! the line-of-sight integration involves the functions jg, ji and j'l. As shown in Ref. 
however, the dependence on j' £ and j" can be eliminated by integrating by part. (We have corrected 
the obvious typos in eq. (12b) in Ref. 

Since is evaluated numerically and it typically oscillates quickly, its sampling is time consum- 
ing. Indeed, in modern fast CMB codes - such as CAMB, CLASS, CMBEASY - the integral ([32]) 
is computed by sampling A^ using a step size in k that can be typically much larger than the 
oscillation period in A^. For instance, in Fig. [2] we show an example of A^ for a fixed t = 300. 
A typical sampling scheme is shown by the red solid triangles in the upper-right panel, which 
zooms-in part of the figure. According to Parseval's theorem, if 'P s (ln/c) is a smooth function, such 
sparse sampling of A^ is enough. 

However, when T , s (k) has local sharp features, the minimum sampling distance should be de- 
termined by the larger between the typical relative width (the width measured in Ink) of the 
oscillations in A^ and 5 Ink, the typical relative width of the features in V s {k). The former is 
about 10 -4 — 10 -3 , while the latter is model-dependent. For instance, if our goal is to sample 
features with width <51nA; ~ 10 -3 , the required sampling frequency is typically ~ 20 to 100 times 
higher than that used for a smooth V s (k). Furthermore, as we wish to compute the Cg's for each £ 
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rather than interpolating it over few tens of £'s, the total time consumption will be again multiplied 
by a factor of ~ 10 — 50. The naively estimated total time consumption is hence ~ 10 3 times more 
than that in the smooth- V s (k) case. A final complication is due to the fact that, if all the transfer 
functions and the precomputed jg(x) tables are to be stored, one has also to face a memory barrier 
that cannot be easily bypassed. For these reasons, simply increasing the I and k resolution in 
standard codes such as CAMB, CLASS or CMBEASY, will not be efficient enough for the purpose 
of scanning the whole parameter space. 

The algorithm can be significantly improved, however, if we notice that the output S(k, r) from 
the Boltzmann code is a 2D matrix in k-r space. If jg [k(ro — r)] is also a precomputed 2D matrix 
with the same structure, the integration (|33p can be obtained by taking the inner product of the 
two matrices. Modern Fortran90 compilers can optimize such operation and make the computation 
much faster. The difficulty, however, is that the jg [k(ro — r)] matrices for all ts will occupy too 
much memory (can be up to a few tens of Giga bytes in the worst scenario). Our solution is then 
to only store the matrices for two neighboring £'s and update them using the recurrence relation 
of spherical Bessel functions. 

Let us describe our strategy. We first compute two neighboring Cg's by brute force. Two 
matrices of spherical Bessel functions jg+\[k{TQ — r)] and jg[k(ro — r)] are stored in the memory 
for each (k, r) indices. Then we compute Cg-\. To do that, we update the jg + \ matrix to the jt-\ 
matrix using the recurrence relation 

2£ + 1 

je-i(x) = jg{x) -jg + i(x) . (34) 

x 

Again, using jg and jg_\ we then calculate jg-i and hence Cg^2- This downward iteration is very 
stable for a few tens of steps, after which we need to recompute another couple of neighboring Cg's 
and iterate downward again. 

The initial neighboring jg's are calculated using precomputed 25-th order Chebyshev fitting 
formulas. (For the rapidly oscillating part at i > I, the phase and amplitude of oscillations are 
fitted using Chebyshev polynomials.) Chebyshev fitting is slightly slower than the cubic-spline 
fitting used in other publicly available CMB codes, but it is more memory-efficient and more 
accurate - it has an accuracy of ~ 10~ 8 - and allows more downward iterative steps. Finally, note 
that the algorithm proposed here is more efficient both CPU-wise and memory-wise, enhancing 
the speed of £-by-i computation of C/s by a factor of ~ 10 to 30. 

For CMB lensing we use the power spectrum approach as described in Refs. 
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B. Testing the Code 

The trivial comparison between CosmoLib and CAMB for smooth- V s models can be found in 



the online documentation at http://www.cita.utoronto.ca/~zqhuang/CosmoLib Here we focus 



on the enhanced CMB integrator that does not assume the smoothness in V s (k). Since this feature 
is not available in other CMB codes, direct numerical comparison is not possible when there is very 
sharp features in V s . Thus, we need to study a model in which we have some theoretical insights. 
An ideal candidate is the axion monodromy inflation model, where the primordial power spectrum 
disp lays sinusoidal oscillations superimposed to a smooth power spectrum. It can be written as 



38] 



= a Att 



r ( Hk/K) 
1 + dn s cos I £Ynk ^ 



(35) 



where A s and n s are the amplitude and tilt, respectively. The parameter 5\nk describes the width 
of the oscillations in V s (k), while 5n s gives their relative amplitude. The pivot scale k* is chosen 
to be 0.05Mpc _1 in our computation. 

We compute the CMB temperature power spectrum using the enhanced CMB integrator in 
CosmoLib and compare the results to the smooth-T , s (k) case. The relative difference between the 
non-smooth (for a series of 5 Ink) and the smooth model is shown in Figure [3l For 6 In k = 0.1 and 
5 In k = 0.03 we compare the results to CAMB output (both with lensing) and find good agreement. 
The CAMB outputs are obtained by a straightforward modification of CAMB, i.e., increasing the 
£ sampling frequency in the input file and increasing the k sampling frequency in the source code. 
For 5\nk < 0.01 the simple modification of CAMB fails due to insufficient memory to store the 
transfer functions. 

For 5 In k <C 1, the amplitude of oscillations in the CMB angular power spectrum (right-hand 
panels) is smaller than that in V s {k) (left-hand panels). This suppression is generic when a 3D 
power spectrum is projected to a 2D one, even though in the CMB case it is further complicated 



by the finite duration of recombination and the recombination physics |74| . As shown in 58(] , when 
the frequency of oscillations is constant in Ink, such as in eq. (I35p . the relative suppression is given 
by ~ V^ln/c, as confirmed by the examples shown in Figure El Moreover, for 5 Ink < 0.01, in 
addition to the projection effect, CMB lensing also significantly smears out the oscillations in Cp at 
high I > 2000. While for 5\nk = 0.1, the lensing smearing effect is almost negligible. See [74I. u\ 
for more detailed discussions about the lensing smearing effect. Finally, note that, although the 
oscillations in Ci are damped, they maintain the same relative width of those of the left-hand 
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panels, i.e., 5ln£ = Sink where £ > (Sink) -1 . At low I where I < 1/5 Ink the oscillations in k 
space disappear in t space due to the discreteness of 1. 

This discussion shows that the enhanced CMB integrator can accurately compute the oscillations 
in CMB to ACg/Cg < 10 -3 . This does not mean, however, that the total Cg is accurate to 10 -3 . 
The Cg power spectrum can be systematically biased at subpercent level due to, e.g., recombination 



uncertainties [76[] . Understanding and eliminating these theoretical errors is important if we want 
to extract generic features in V s {k) with 10 -3 accuracy. On the other hand, if we are only interested 
in a model predicting a specific feature in Cg that cannot be mimicked by other effects, we can 
focus only on the relative difference in Cg. 

IV. THE FORECAST TECHNIQUES 

CosmoLib uses Fisher matrix analysis and MCMC method to forecast the constraints on cos- 
mological parameters for future CMB, LSS and SN experiments. In this section we discuss the 
modeling of the likelihoods and the parameter estimation methods. 

A. The likelihoods 

1. CMB simulation 



Given a likelihood function 



X 2 can be approximated by 



77 



r w e define x 2 = — 21n£. For a nearly full-sky CMB experiment 



78] 



X 2 = J2 (2^+U/sky 



pBB / nBB \ pTTpEE , pEEpTT npTEpTE I pTTpEE _ (pTE\2 

C BB + I qBB I + C JT C EE _ {C TE ? + ln \fiTT£EE _ 



TE\2 



(36) 

where ^ m ; n and £ mSLX are suitable cutoffs that are determined by the observed fraction of sky / s k y 
and the beam resolution of the experiment. In this formula, Cf Y are the model-dependent the- 
oretical angular power spectra (including the noise contributions) for the temperature, E and B 
polarizations and their cross-correlations, with X, Y = {T, E, B}. We compute the noise contribu- 
tion Ng assuming Gaussian beams. The mock data C¥ Y are Cf Y for the fiducial model. 



We use the model introduced in [77j] (and later updated in |78j|) to propagate the effect of 
polarization foreground residuals into the estimated uncertainties on the cosmological parameters. 
For simplicity, in our simulation we consider only the dominant components in the frequency bands 
that we are using, i.e., the synchrotron and dust signals. We assume that foreground subtraction 



15 



can be done correctly down to a level of 5%. (This parameter is adjustable by the user.) 

2. SN simulation 

For the SN simulation, we use the model given by the Dark Energy Task Force (DETF) forecast 
[79! - In this case 

* 2 -E(^y. <37) 

with i going over the SN samples. More specifically, here rrii and rhi are the theoretical expectation 
and observed magnitude of the z-th supernova, respectively. The uncertainty Srrii is computed by 
quadratically adding a peculiar velocity (a user-defined constant) to the intrinsic uncertainty in 
the supernova absolute magnitude (another user-specified constant). 
The apparent magnitude of SN is modeled as 

m = M-fi L z- fi Q z 2 + 5 log 10 + 25 - n S S Qeax . (38) 

The first three terms model the redshift evolution of the absolute magnitude of the supernova peak 
luminosity. In particular, M is a free parameter with a flat prior over —00 < M < +00; for \ir 
and Gaussian priors are applied. The widths of the Gaussian priors are user-input parameters. 
Finally, given that the nearby samples are likely to be a collection from many other experiments, 
an offset — /x 5 <5 ncar , where <5 nea r is unity for the nearby samples [z < z nca , r ) and zero otherwise, is 
added to account for the systematics. For fi s we also apply a Gaussian prior with a user-specified 
width. The threshold redshift z ncai is also user-defined. In conclusion, in this model there are four 
nuisance parameters (M, fi , fjS, /j, s ), which we marginalized analytically. 

3. LSS likelihood 



We model the galaxy power spectrum in redshift space as (e.g., 8CH82I] ) 



P g (k,wz) = {b + f f i 2 ) 2 D 2 (z)P m (k)exp(-k 2 f i 2 a 2 r ) , (39) 

where /i is the cosine of the angle between the wave vector k and the line of sight, D{z) is the linear 
growth factor, / = dlnD/d\na is the linear growth rate, P m (k) is the matter power spectrum today 
(at z = 0) and a r parameterizes the effect of small scales velocity dispersion and redshift errors as 
explained below. The matter power spectrum P m (k) is computed using Poisson's equation, that 
is, P m (k) = 4k^ k \ 2 /(9H 4 n 2 J. 
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The term f/j? accounts for the redshift distortions due to the large-scale peculiar velocity field 



80l | . which is correlated with the matter density field. The exponential factor on the right-hand 



side accounts for the radial smearing due to the redshift distortions that are uncorrelated with the 
LSS. In particular, we consider two contributions. The first is due to the redshift uncertainty of the 
spectroscopic galaxy samples which is estimated to be a z = 0.001(1 + z) [83]. (In CosmoLib the 
user is allowed to change this value.) The second comes from the Doppler shift due to the virialized 
motion of galaxies within clusters, which typically has a pairwise velocity dispersion a g of the order 
of few hundred km/s. This can be parameterized as ^|(1 + z ) |£2|]- The two contributions are 
quadratically added together 

^ 2 = ^f W + °2 /2 ) ' (40) 

where H{z) is the Hubble parameter. 

Practically, neither the redshift measurement nor the virialized motion of galaxies can be pre- 
cisely modeled. In particular, the radial smearing due to peculiar velocity is not necessarily close 
to Gaussian. Thus, eq. (|39p should not be used for wavenumbers k > - > where the radial 
smearing effect is important. We introduce a UV cutoff fc max as the smallest value between - 
and t^t, where R is chosen such that the r.m.s. linear density fluctuation of the matter field in a 
sphere with radius R is 0.5. 

The survey volume is split into n z redshift bins from z m ; n to z max , with all these parameters to 
be specified by the user. The number density of galaxies that can be used is n = en b s , where e is 
the fraction of galaxies with measured redshift to be specified by the user. Due to the high accuracy 
of the spectroscopic redshift and the width of the bins, we ignore the bin-to-bin correlations and 



write x 2 as 



2 _ I -fg.modcl -Pg,fiducial \ 

X L> V AP a . fiducial • { ' 

k,fi,z bins a ' 

As on large scales the matter density field has, to a very good approximation, Gaussian statistics 
and uncorrelated Fourier modes, the band-power uncertainty is given by 



2 



84] 



AP g 



2(2vr) 



3 iV2 



(2Trk 2 dkdfi) (4irr 2 f sky dr) 

where r is the comoving distance given, for a flat FRW universe, by r(z) = f cdz' /H(z'). The 
second term in the parenthesis is due to shot noise, under the assumption that the positions of the 
observed galaxies are generated by a random Poisson point process. In practice n is not known 
a priori and is calibrated by galaxies themselves. The imperfect knowledge of n can bias P g on 
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the scale of the survey 84]. This is taken into account by using an IR cutoff k m - m ~ Gpc -1 . This 
is chosen such that k® in = 2ir/V^ 3 , where V{ is the comoving volume of the i-th. (i = 1,. . . ,n z ) 
redshift slice. Finally, the user has to specify the binning scheme for k and fi. For k we allow 
uniform binning in In A; or in k. For ji only uniform binning in \x is allowed. 

In the special case where P m (k) has sharp features, we must consider the smearing effect due to 
the fact that we are only observing a finite volume. This effect is approximated by replacing P m (k) 
in (|39p with its convolution with a Gaussian window, where the width of the Gaussian window aw 
has been chosen to be 

v / 21n2 /4tt\ 1/3 . , 

ow = ( y J fcmin ^ 0.302 k min . (43) 

In such a way, the real-space representation of the window, if cut off at its half-height, contains 
the same volume as that of the redshift bin. The fact that % is smaller than k m { n allows us to 
neglect the overlap between window functions centered around neighboring values of k. 

B. Parameter Estimation 

1. Fisher Matrix Analysis 

In general, the likelihood can be written as 

ln£(p;p fid ) = ~ ]d(p) - d(p fid )] r C r - 1 (p; Pfid ) ]d(p) - d(p fid )] , (44) 

where d is the data vector, pfid the fiducial parameter vector, p the parameter vector for which 

one wants to evaluate the likelihood, and C~ 1 (p;pfi ( j) the covariance matrix. 

The fisher matrix for pi, pj (two components of p) is then 

d 2 lnC dd(p) x 0d(p) 

= — C (p;pm)^ — 



13 dpidpj 



d Pi dp 3 



(45) 
p=Pfld 



where the partial derivatives can be evaluated numerically: 

P = [d(pi,P2,.. ■ ,Pi + Apt, ■ ■ ■ ,p n ) - d(pi,p 2) . • • ,Pi - Api, . . . ,p n )] ■ (46) 

The stepsize Api is initially supplied by the user, and then adjusted by the software in such a 
way that the variation of x 2 is of 0(1) when pi is varied by Apj. By doing this, we have assumed 
that the likelihood is approximately Gaussian in the proximity of pf^ where the variation of x 2 
is < O(l). If the likelihood is highly non-Gaussian, Fisher matrix analysis does not give reliable 
estimations of the error bars of parameters. In this case, one should use the MCMC method to 
fully explore the structure of the likelihood. 
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2. MCMC method 

CosmoLib has an independent MCMC engine using the Metropolis-Hastings algorithm. The 
traditional approach is to define the proposal density Q(x;x') (the probability of walking from x 
to x' in the parameter space) using a roughly estimated covariance matrix C e 



Q(x; x') oc exp 



■^(x-x')V(x-x') 



(47) 



Convergence can be achieved quickly if C e is close the posterior covariance matrix of x. 

However, sometimes we need to treat models where the likelihood periodically depends on some 
phase parameters. Here we take the axion monodromy inflation model for example. The likelihood 
L is a periodic function of the axion phase parameter ip, 

£(P,p)=£(P,P + 2tt) , (48) 

where we have used P to represent the collection of other parameters. If ip is not well constrained, 
we will obtain multi-branches in the posterior, i.e., for a fixed value of p and a chosen confidence 
level, the allowed values of P locate in well separated regions in the parameter space. 

Intuitively the separated regions can be more efficiently explored by restricting the range of p 
to one period and proposing with wrap-around or, in a more rigorous language, by using a periodic 
proposal density. For x = (P, ip) and x' = (P' , ip'), we use 



Q(p,<p;P',<p') oc ^2 ex P 



" 2 ( X ~~ X n) T Ce 1 ( x — x n 



(49) 



where x' n = (P', p' -\-2nir). The estimation of covariance matrix, C e , is practically computed with a 
trial run that is terminated before the multi-branches of the posterior are explored by the random 
walk. 

We find that the periodic proposal density (|49|) leads to significant improvement of the conver- 
gence. For the axion monodromy model, it takes about 5-10 times longer to achieve convergence 
using (|37|) than using (|4"9"j). 

The output chains in CosmoLib have the same format as those in CosmoMC 



. The chains can 



hence be directly analyzed using the GetDist tool in CosmoMC. For completeness, an independent 
postprocessing tool is supplied in CosmoLib to analyze and visualize the marginalized posterior of 
parameters. In the online documentation the reader can find the instructions on how to use this 
tool. 
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3. Oscillations in the Current CMB Data? 



Recently a hint of the axion monodromy cosine oscillations (see eq. (I35j) ) in WMAP-7yr 10l.l8E| 



and ACT CMB data 



87[ has been claimed in Ref . 56J] . Ref. [57( confirms the finding that x can be 



significantly improved in some regions of parameter space where oscillations in the primordial power 
spectrum are assumed. In this section we use CosmoLib to constrain the axion monodromy model 
with the same data sets. We find that when the CMB power spectrum is accurately computed 
and rigorous statistical method is used, there is no detectable axion monodromy oscillations in the 
CMB data. 



In Refs. |56l ] the authors used their modified CAMB to compute the CMB power spectrum. As 
discussed in previous sections, such a modification is not trivial for <51n k < 10 -2 . Since the best-fit 
5 In k found in Ref. 



51 in Ref. 



561 ] and equation 



561 ] is small - 5 In k ~ 0.005 (derived from Table III of Ref. 
58l ] ) , it is necessary to exam the numerical accuracy of the modified CAMB used in 
For Sink ~ 0.005, the modulation period in In A; is T\ n k = 2ir5lnk ~ 0.03. In the CMB power 
spectrum one should see same modulation period in ln£, i.e., T\ n £ = T\ n ^ w 0.03. Thus, from 
£ = 1000 to £ = 1200 there should be about 7 oscillations in Cg. However, in Figure 5 of Ref. 
the number of oscillations in Cp between £ = 1000 and £ = 1200 are much more than 7. This implies 



56] 



that the "modulations" in Cg shown in Ref. 



561 ] may just be numerical noises. In Figure H] we show 



the CMB temperature angular power spectrum computed with CosmoLib, where the parameters 
are chosen to be close to the ones used in Figure 5 of Ref. [56J. Qualitative difference can be seen 



between the two figures. The Cg spectrum computed using CosmoLib presents clear modulations 
that agrees with the 5 In k value, while the modified CAMB used in Ref. 
expected modulations. 



561 ] failed to produce the 



In Ref. 58|] we pointed out that, a significant improvement of x 2 does not necessarily imply a 
detection of models with periodic oscillations, which typically has a spiky likelihood that is highly 
non-Gaussian. A rigorous treatment is to compute the marginalized probability of the amplitude of 
oscillations 5n s . The marginalization should be done in such a way that all the other cosmological 
and nuisance parameters are allowed to vary. A detection of monodromy oscillations should not 
be claimed unless 5n s = is excluded by the data. We did the full marginalization using MCMC 
method. The CMB power spectra are computed using the accurate integrator in CosmoLib. The 
marginalized 68.3% and 95.4% confidence level posterior contours are shown in Figure [5j 
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V. CONCLUSIONS 



We introduced the numerical package CosmoLib and focused on its features that are comple- 
mentary to other numerical codes. The major advantage of CosmoLib is that it can accurately 
compute CMB angular power spectrum for inflationary models that predict sharp features in the 
primordial power spectrum of metric perturbations. This is not available in any other publicly 
available CMB codes. 

CosmoLib can calculate the relative fluctuations in Cg to accuracy ~ 10~ 3 . Because of cosmic 
variance, we cannot measure Cg to this accuracy if all Cg are treated independently. However, 
our purpose is to use CosmoLib to study specific models where the degrees of freedom in the Cg 
spectrum is small. In other words, if we assume a specific model (such as the axion monodromy 
model), the relative error in Cg can be constrained to a level that is well below cosmic variance. 
In the naive limit where the Cg spectrum is controlled by a single scaling parameter s, that is, 
Cg = sC^fiduciab we can constrain Cg to a relative accuracy 1/VN ~ 1/Anax, where iV = ^2g(2£ + 
1) ~ ^max i s the total number of measured spherical harmonic modes. For a future experiment that 
measures Cg to cosmic variance for £ up to a few thousands 88(, the aforementioned 10 -3 relative 
accuracy is necessary. 

While a straightforward (but not optimized) modification of CAMB and CLASS to use non- 
smooth P(k) seems to be trivial, in practice it is often limited by the available memory and tolerable 



561 ] produces numerical noises 



computation time. We pointed out that the modified CAMB in Ref. 
instead of the expected modulation in Cg spectrum. Repeating the computation in Ref. [! 
CosmoLib and the same data sets (WMAP + ACT), we found no detection or hint of axion 
monodromy model in the current CMB data. 

This forecast toolkit contains a fisher matrix calculator, a MCMC engine, a postprocessing tool 
for chain analysis, and likelihoods for future CMB, galaxy survey, and supernova observations. 
The MCMC engine has an option of using a periodic proposal density, which can significantly 
accelerate the convergence of the chains in the case where the likelihood is a periodic function of 
some parameters. Although the likelihood models in CosmoLib are likely to be too simple for real 
experiments with complicated specifications, they provide a quick estimation of the performance 
of future CMB/LSS/SN experiments, for which the details of specifications are not yet known. We 
are planning to include more likelihoods for, e.g., weak lensing experiments in future releases. 
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FIG. 3. The differences in \n.V s (left panels) or \nCj T (right panels) between a fiducial axion monodromy 
model with In (l0 lo A s ) = 3.027, n s = 0.975, amplitude of cosine modulation Sn s = 0.01, phase tp = and 
a smooth power-law spectrum with the same A s and n s . For the top to bottom a series of 8 In k = 0.1, 
0.03, 0.01 are used, respectively. The r rocom b in the x-axis legend of left panels is the conformal time at 
recombination (z ps 1100). For Sink = 0.1 and 0.03 the results are compared to CAMB outputs. For 
Sink = 0.01 a simple modification of CAMB cannot be applied as too much memory is required to store 
the transfer functions for all (£, k) pairs. 



26 




FIG. 4. The CMB angular power spectrum for axion monodromy model with S In k = 0.005, Sn s = 0.18. The 
other cosmological parameters are tt b0 h 2 = 0.0223, Q c0 h 2 = 0.1119, 6 = 1.041, r re = 0.0884, n s = 0.975, 
ln(10 10 A s ) = 3.04. Modulation of Cg is uniform in ln£ and is almost invisible at high-€ due to lensing 
smearing. This should be compared to Fig. 5 in Ref. |56j, where the random fluctuations in Ci implies 



insufficient numerical accuracy of the modified CAMB used therein. 




5 Ink 

FIG. 5. The marginalized 68.3% and 95.4% confidence level contours of 8n s and Sink for axion monodromy 
model. WMAP-7yr and ACT data are used. CMB angular power spectrum are computed up to I = 4000 
with CMB lensing effect included. Uniform priors 0.003 < 5 In k < 0.2 and < n s < 0.2 are used. 
No detection of axion monodromy oscillations are found since zero amplitude of oscillations (5n s = 0) is 
consistent with the data. 



